# data wrangling
library(tidyverse, warn.conflicts = FALSE, quietly = TRUE)
library(tidyr, warn.conflicts = FALSE, quietly = TRUE)
library(dplyr, warn.conflicts = FALSE, quietly = TRUE)
library(ggplot2)
library(ggpubr)
library(sjPlot)
library(kableExtra)
library(xtable)
library(plotly)
library(ggthemes)
library(ggbeeswarm)
# linear model
library(lme4)
library(afex) # mixed effect model
library(permuco)
library(remotes)
library(permutes)
library(coin)
rm(list=ls())
set.seed(0)
Load Maximum LogLikelihood data
d <- read.csv("./max_loglikelihood.csv", sep = ",") %>%
select(HCPID, best_model, maxLL_diff,
Punishment_MostlyPunishment.subj,
Reward_MostlyPunishment.subj,
Reward_MostlyReward.subj,
Punishment_MostlyReward.subj) %>%
pivot_longer(Punishment_MostlyPunishment.subj:Punishment_MostlyReward.subj,
names_to = "Condition",
values_to = "Pswitch") %>%
mutate(
block_type = case_when(
Condition == "Reward_MostlyPunishment.subj" ~ "Loss",
Condition == "Punishment_MostlyPunishment.subj" ~ "Loss",
Condition == "Reward_MostlyReward.subj" ~ "Reward",
Condition == "Punishment_MostlyReward.subj" ~ "Reward",
TRUE ~ "NA"),
trial_type = case_when(
Condition == "Punishment_MostlyReward.subj" ~ "Loss",
Condition == "Punishment_MostlyPunishment.subj" ~ "Loss",
Condition == "Reward_MostlyReward.subj" ~ "Reward",
Condition == "Reward_MostlyPunishment.subj" ~ "Reward",
TRUE ~ "NA"),
)
#d.trial = d %>%
# group_by(HCPID, trial_type, best_model) %>%
# summarise(Pswitch = mean(Pswitch))
We could fit data using glmer
trial_type, block_typebest_modelHCPIDFitting data to mixed effect linear model, best_model label is significant
summary(m <- glmer(data=d, Pswitch ~ d$best_model + (trial_type+block_type|HCPID),
family = "binomial"))
Let’s try other packages.. apex library
Significant fixed effects:
best_model
trial_type
2 way interaction effects:
best_model : block_type
best_model : trial_type
3 way interaction effects:
best_model:block_type:trial_type#apex library
mixed_anova <- aov_ez(
id = "HCPID",
dv = "Pswitch",
data = d,
between = "best_model",
within = c("block_type", "trial_type")
)
print(mixed_anova)
## Anova Table (Type 3 tests)
##
## Response: Pswitch
## Effect df MSE F ges p.value
## 1 best_model 1, 197 0.10 43.55 *** .099 <.001
## 2 block_type 1, 197 0.02 1.22 <.001 .270
## 3 best_model:block_type 1, 197 0.02 6.92 ** .004 .009
## 4 trial_type 1, 197 0.05 7.33 ** .009 .007
## 5 best_model:trial_type 1, 197 0.05 41.10 *** .046 <.001
## 6 block_type:trial_type 1, 197 0.03 0.43 <.001 .515
## 7 best_model:block_type:trial_type 1, 197 0.03 12.16 *** .009 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
ggplotly(afex_plot(mixed_anova, x = "best_model", trace = c("trial_type"), #error = "none",
mapping = c("color", "shape"), data_arg = list(color = "gray", shape = 1),
line_arg = list(color = c('tomato2', 'tomato2', 'forestgreen', 'forestgreen'), size=2),
data_plot = T, point_arg = list(size = 5,
color = rep(c('tomato2', 'forestgreen'), each=2),
shape = rep(c('x', '.'), each=2))) +
labs(title = "Behaviror: Trial effect between two models",
x = "ACT-R Model Type", y = "Probability of Switch") +
theme_bw())
ggplotly(afex_plot(mixed_anova, x = "best_model", trace = c("block_type"), #error = "none",
mapping = c("color", "shape"), data_arg = list(color = "gray", shape = 1),
line_arg = list(color = c('tomato2', 'tomato2', 'forestgreen', 'forestgreen'), size=2),
data_plot = T, point_arg = list(size = 5,
color = rep(c('tomato2', 'forestgreen'), each=2),
shape = rep(c('x', '.'), each=2))) +
labs(title = "Behaviror: Block effect between two models",
x = "ACT-R Model Type", y = "Probability of Switch") +
theme_bw())
Similarly, fit data into mixed effect model
m1 <- mixed(Pswitch ~ best_model + (trial_type + block_type | HCPID), d,
method = "LRT",
all_fit = TRUE,
family = binomial)
nice(m1)
## Mixed Model Anova Table (Type 3 tests, LRT-method)
##
## Model: Pswitch ~ best_model + (trial_type + block_type | HCPID)
## Data: d
## Df full model: 8
## Effect df Chisq p.value
## 1 best_model 1 63.59 *** <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
xtable(anova(m1))
## % latex table generated in R 4.0.3 by xtable 1.8-4 package
## % Fri May 13 22:04:54 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{lrrrr}
## \hline
## & Df & Chisq & Chi Df & Pr($>$Chisq) \\
## \hline
## best\_model & 7 & 63.59 & 1 & 0.0000 \\
## \hline
## \end{tabular}
## \end{table}
summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Pswitch ~ best_model + (trial_type + block_type | HCPID)
## Data: data
## Control: ctrl
##
## AIC BIC logLik deviance df.resid
## 1027.5 1064.9 -505.7 1011.5 788
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.20028 -0.22055 0.05622 0.42645 1.52324
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## HCPID (Intercept) 0.000e+00 0.000e+00
## trial_type1 5.095e-15 7.138e-08 NaN
## block_type1 1.981e-14 1.407e-07 NaN -0.11
## Number of obs: 796, groups: HCPID, 199
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.23828 0.07699 -3.095 0.00197 **
## best_model1 -0.60340 0.07699 -7.837 4.6e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## best_model1 -0.211
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Next, we could use apex package to perform permutation test.
d.block = d %>%
group_by(HCPID, best_model, block_type) %>%
summarise(Pswitch = mean(Pswitch), .groups = "drop") %>%
spread(key = block_type, value = Pswitch) %>%
mutate(block_effect = Reward - Loss) %>%
select(HCPID, best_model, block_effect) %>%
spread(key = best_model, value = block_effect)
d.trial = d %>%
group_by(HCPID, best_model, trial_type) %>%
summarise(Pswitch = mean(Pswitch), .groups = "drop") %>%
spread(key = trial_type, value = Pswitch) %>%
mutate(trial_type = Reward - Loss) %>%
select(HCPID, best_model, trial_type) %>%
spread(key = best_model, value = trial_type)
d.block.m1 <- d.block %>% select(m1) %>% drop_na()
d.block.m2 <- d.block %>% select(m2) %>% drop_na()
d.trial.m1 <- d.trial %>% select(m1) %>% drop_na()
d.trial.m2 <- d.trial %>% select(m2) %>% drop_na()
compare.2.vectors(d.block.m1$m1, d.block.m2$m2, paired = FALSE, na.rm = FALSE,
tests = c("parametric", "nonparametric"), coin = TRUE,
alternative = "two.sided",
wilcox.exact = NULL, wilcox.correct = TRUE)
## $parametric
## test test.statistic test.value test.df p
## 1 t t -2.630521 197.0000 0.009199408
## 2 Welch t -2.661803 152.8869 0.008603209
##
## $nonparametric
## test test.statistic test.value test.df p
## 1 stats::Wilcoxon W 3604.000000 NA 0.01318405
## 2 permutation Z -2.592059 NA 0.00959000
## 3 coin::Wilcoxon Z -2.480040 NA 0.01244000
## 4 median Z -2.407694 NA 0.01776000
compare.2.vectors(d.trial.m1$m1, d.trial.m2$m2, paired = FALSE, na.rm = FALSE,
tests = c("parametric", "nonparametric"), coin = TRUE,
alternative = "two.sided",
wilcox.exact = NULL, wilcox.correct = TRUE)
## $parametric
## test test.statistic test.value test.df p
## 1 t t 6.411169 197.0000 1.045368e-09
## 2 Welch t 6.534571 156.1739 8.553576e-10
##
## $nonparametric
## test test.statistic test.value test.df p
## 1 stats::Wilcoxon W 6826.500000 NA 7.72669e-09
## 2 permutation Z 5.846384 NA 0.00000e+00
## 3 coin::Wilcoxon Z 5.775596 NA 0.00000e+00
## 4 median Z 4.338846 NA 2.00000e-05
perm.lmer(Pswitch ~ best_model + (trial_type + block_type | HCPID), data=d, family = binomial(), type = "regression")
## Factor LRT beta z p
## 1 (Intercept) 187.3028 -0.4147827 -4.575609 NA
## 2 best_modelm2 115.7080 0.6174496 4.139432 0
We could also simulate mean difference and calculate the p-value by creating our own permutation test
Block effect: Pswitch(Reward Block) - Pswitch(Loss Block)
Trial effect: Pswitch(Reward Trial) - Pswitch(Loss Trial)
Mean Difference of Block Type: Block effect(M1) - Block effect(M2)
Mean Difference of Trial Type: Trial effect(M1) - Trial effect(M2)
shuffle_labels <- function(d) {
res <- d %>%
select(HCPID, best_model) %>%
unique()
res$best_model <- sample(res$best_model, replace = T)
res <- d %>% select(-best_model) %>%
left_join(res, by = "HCPID")
return (res)
}
block_effect <- function(dat) {
res = dat %>% group_by(best_model, block_type) %>%
summarise(Pswitch = mean(Pswitch), .groups = 'drop') %>%
pivot_wider(names_from = block_type, values_from = Pswitch) %>%
mutate(block_effect.RL = Reward-Loss) %>%
select(best_model, block_effect.RL) %>%
pivot_wider(names_from = best_model, values_from = block_effect.RL) %>%
mutate(diff.m1m2 = m1-m2)
return (res)
}
trial_effect <- function(d) {
res = d %>% group_by(best_model, trial_type) %>%
summarise(Pswitch = mean(Pswitch), .groups = 'drop') %>%
pivot_wider(names_from = trial_type, values_from = Pswitch) %>%
mutate(trial_effect.RL = Reward-Loss) %>%
select(best_model, trial_effect.RL) %>%
pivot_wider(names_from = best_model, values_from = trial_effect.RL) %>%
mutate(diff.m1m2 = m1-m2)
return (res)
}
my_block_effect <- block_effect(d)
my_trial_effect <- trial_effect(d)
# simulate random block effect
n <- 5000
simulated_block_effect = c()
simulated_trial_effect = c()
for (x in 1:n) {
d0 <- shuffle_labels(d)
curr.block <- block_effect(d0)
curr.trial <- trial_effect(d0)
simulated_block_effect <- c(simulated_block_effect, curr.block$diff.m1m2)
simulated_trial_effect <- c(simulated_trial_effect, curr.trial$diff.m1m2)
}
Calculate one-side and two-side p-value
# Let's compute the p-value
if (my_block_effect$diff.m1m2 < 0) {
pvalue.block <- sum(my_block_effect$diff.m1m2 > simulated_block_effect) / n
} else {
pvalue.block <- sum(my_block_effect$diff.m1m2 <= simulated_block_effect) / n
}
paste("Block Type: The one-sided p-value is ", round(pvalue.block,3))
## [1] "Block Type: The one-sided p-value is 0.005"
paste("Block Type: The two-sided p-value is ", 2 * round(pvalue.block,3))
## [1] "Block Type: The two-sided p-value is 0.01"
if (my_trial_effect$diff.m1m2 < 0) {
pvalue.trial <- sum(my_trial_effect$diff.m1m2 > simulated_trial_effect) / n
} else {
pvalue.trial <- sum(my_trial_effect$diff.m1m2 <= simulated_trial_effect) / n
}
paste("Trial Type: The one-sided p-value is ", round(pvalue.trial,3))
## [1] "Trial Type: The one-sided p-value is 0"
paste("Trial Type: The two-sided p-value is ", 2 * round(pvalue.trial,3))
## [1] "Trial Type: The two-sided p-value is 0"
#pvalue.greater <- sum(my_block_effect$diff.m1m2 > simulated_block_effect) / n
#pvalue.smaller <- sum(my_block_effect$diff.m1m2 <= simulated_block_effect) / n
#paste("One-sided p-value (greater) is", pvalue.greater)
#paste("One-sided p-value (smaller) is", pvalue.smaller)
Plot the distribution of mean difference of Block Type between declarative model vs. procedural model
gghistogram(simulated_block_effect, add = "mean", add_density = F, bins = 50, rug = T, fill = 'gray',
title = "Histogram of simulated mean difference: Block Type", xlab = "mean_difference", ylab="simulated counts") +
geom_vline(aes(xintercept=my_block_effect$diff.m1m2), linetype="dashed", color = "red") +
annotate(geom = "text", x = my_block_effect$diff.m1m2*0.8, y = 200,
label = paste("observed p = ", 2 * round(pvalue.block,3)))
Plot the distribution of mean difference of Trial Type between declarative model vs. procedural model
gghistogram(simulated_trial_effect, add = "mean", add_density = F, bins = 50, rug = T, fill = 'gray',
title = "Histogram of simulated mean difference: Trial Type", xlab = "mean_difference", ylab="simulated counts") +
geom_vline(aes(xintercept=my_trial_effect$diff.m1m2), linetype="dashed", color = "red") +
annotate(geom = "text", x = my_trial_effect$diff.m1m2*0.8, y = 200,
label = paste("observed p = ", 2 * round(pvalue.trial,3)))
There are different ways of calculating maxLL. Next we could explore which way gives us most distinguishable model labels
Group by Block Type and Trial Type
Group by Block Type only
Group by Trial Type only
d.block <- read.csv("./max_loglikelihood_block.csv", sep = ",") %>%
select(HCPID, best_model, maxLL_diff,
MostlyPunishment.subj,
MostlyReward.subj) %>%
pivot_longer(MostlyPunishment.subj:MostlyPunishment.subj,
names_to = "Condition",
values_to = "Pswitch") %>%
mutate(
block_type = case_when(
Condition == "MostlyReward.subj" ~ "Reward",
Condition == "MostlyPunishment.subj" ~ "Loss",
TRUE ~ "NA"))
d.trial <- read.csv("./max_loglikelihood_trial.csv", sep = ",") %>%
select(HCPID, best_model, maxLL_diff,
Reward.subj,
Punishment.subj) %>%
pivot_longer(Reward.subj:Punishment.subj,
names_to = "Condition",
values_to = "Pswitch") %>%
mutate(
block_type = case_when(
Condition == "Reward.subj" ~ "Reward",
Condition == "Punishment.subj" ~ "Loss",
TRUE ~ "NA"),
trial_type = case_when(
Condition == "Reward.subj" ~ "Reward",
Condition == "Punishment.subj" ~ "Loss",
TRUE ~ "NA"))
Below plot shows the differences of MaxLL between two models. Different color indicates different ways of calculating MaxLL.
For Both, the maxLL was calcualted with aggregated Pswitch by both block type and trrial type
For block_type, the maxLL was calcualted with aggregated Pswitch only by block type
For trial_type, the maxLL was calcualted with aggregated Pswitch only by trial type
x-axis is the absolute difference between maxLL of model1 (declarative model) and maxLL of model2 (procedural model).
Take-away: If only look at block effect (blue density curve), two model are more distinguisable from each other.
maxLL_diff = data.frame(maxLL_type=c(rep("both", length(d$maxLL_diff)),
rep("block_type", length(d.block$maxLL_diff)),
rep("trial_type", length(d.trial$maxLL_diff))),
maxLL_diff = c(d$maxLL_diff,
d.block$maxLL_diff,
d.trial$maxLL_diff)) %>%
mutate(value_type=factor(maxLL_type, levels = c("both", "block_type", "trial_type")))
ggdensity(data=maxLL_diff, x="maxLL_diff", fill = "maxLL_type", palette = "Set1", add = "mean",
color = "darkgray", alpha = .5,
xlab = "max loglikelihood abs(model1 - model2)",
title = "Density Plot of Max LogLikelihood Differences") +
theme_pander()
Let’s check whether three different ways give us reliable individual strategy selection.
best_models = d.trial %>% select(HCPID, best_model) %>% unique() %>%
left_join(d.block %>% select(HCPID, best_model) %>% unique(), by = "HCPID", suffix = c(".trial", ".block")) %>%
left_join(d %>% select(HCPID, best_model) %>% unique(), by = "HCPID") %>%
mutate(m1_count = if_else(best_model=='m1', 1, 0) +
if_else(best_model.trial=='m1', 1, 0) +
if_else(best_model.block=='m1', 1, 0)) %>%
arrange(m1_count)
Unfortunately, if we use block type or trial type solely to calcualte maxLL, it simply believes all subjects use m1. This is clearly not what we wants.
best_models %>%
count(best_model)
## # A tibble: 2 x 2
## best_model n
## * <chr> <int>
## 1 m1 126
## 2 m2 72
best_models %>%
count(best_model.trial)
## # A tibble: 2 x 2
## best_model.trial n
## * <chr> <int>
## 1 m1 196
## 2 m2 2
best_models %>%
count(best_model.block)
## # A tibble: 1 x 2
## best_model.block n
## * <chr> <int>
## 1 m1 198
d %>% select(HCPID, best_model, maxLL_diff) %>%
arrange(best_model, desc(maxLL_diff))
## # A tibble: 796 x 3
## HCPID best_model maxLL_diff
## <chr> <chr> <dbl>
## 1 164131_fnca m1 7.50
## 2 164131_fnca m1 7.50
## 3 164131_fnca m1 7.50
## 4 164131_fnca m1 7.50
## 5 169141_fnca m1 5.53
## 6 169141_fnca m1 5.53
## 7 169141_fnca m1 5.53
## 8 169141_fnca m1 5.53
## 9 155231_fnca m1 5.16
## 10 155231_fnca m1 5.16
## # … with 786 more rows